{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "406c2750",
   "metadata": {},
   "source": [
    "# Bike Share Rebalancing With Mathematical Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a7e7404",
   "metadata": {},
   "source": [
    "Bike share systems have become an effective commuting method globally for everyday urban dwellers as well as tourists.\n",
    "\n",
    "Citi-Bike in NYC being the largest Bike-Sharing network had 1,588 active stations and 25,575 active bikes in July 2022.\n",
    "\n",
    "Over 3 million rides were completed July 2022 that cover NYC/Hoboken/Jersey City, with around 150,000 active annual members.\n",
    "\n",
    "During rush hours there are many bike stations that have a high demand for bikes, which means their out-flow of bikes is greater than their in-flow in these stations. \n",
    "\n",
    "Meanwhile there are stations that have a high demand for docks (riders return their bikes to these stations) which means their in-flow of bikes is greater than their out-flow.\n",
    "\n",
    "Lack of available bikes or docks in high-demand stations can cause major imbalance in the bike sharing network and result in customer dissatisfaction and lost revenue.\n",
    "\n",
    "To tackle this problem, bikes are relocated between stations to create a balance between supply and demand."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fdb0fdf",
   "metadata": {},
   "source": [
    "## Problem Statement and Solution Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2142854",
   "metadata": {},
   "source": [
    "Using historical Citi-bike data in NYC and Jersey area during July 2022, we like to know:\n",
    "- What is the demand for bikes per hour at each station during the first week of August?\n",
    "- Knowing the demand, how can we minimize loss of sale?\n",
    "\n",
    "Loss of sale is caused by lack of bikes when customers demand them. So, bikes should be transferred from stations with higher in-flow of bikes to those with higher out-flow of bikes. \n",
    "\n",
    "So, first, number of bikes to be added to or removed from each station during each hour should be determined. Then, the physical transfer of bikes between stations should be scheduled. \n",
    "\n",
    "In this notebook, we'll focus on the first part and at the end, discuss how the second part can be solved.\n",
    "We'll use a mixture of Machine Learning (ML) and Mathematical Optimization (MO) to solve this problem. \n",
    "\n",
    "**Solution Approach**\n",
    "The solution approach is comprised of two steps:\n",
    "- **Step 1**: We use the historical Citi-bike data in NYC and Jersey area during July 2022 and use an ML model to predict the number of in-flow and out-flow of bikes per hour at each station for the first week of August. This is done in [predict_bike_flow](predict_bike_flow.ipynb) Notebook.\n",
    "- **Step 2**: We use an MO model to decide how many bikes should be added to or removed from each station during each hour so that the total loss of sale is minimized.\n",
    "\n",
    "To ensure that everyone can run the notebook with the Gurobi restricted license, we reduce the size of the data. To achieve that, we focus on the top 50 stations during the morning rush hours (7 am to 9 am).\n",
    "\n",
    "The top stations are chosen using the PageRank algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed4f7ee5",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Install Required Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af590656",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-08T19:29:24.273068Z",
     "start_time": "2022-11-08T19:29:20.555378Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "%pip install gurobipy\n",
    "%pip install pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8bd6f69",
   "metadata": {},
   "source": [
    "# Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4036b5a7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-16T00:27:52.744267Z",
     "start_time": "2022-11-16T00:27:52.303702Z"
    }
   },
   "outputs": [],
   "source": [
    "import datetime\n",
    "import gurobipy as gp\n",
    "import pandas as pd\n",
    "from gurobipy import GRB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c165d441",
   "metadata": {},
   "source": [
    "# Optimization Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58e8f88d",
   "metadata": {},
   "source": [
    "## Problem Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0816abcb",
   "metadata": {},
   "source": [
    "We want to minimize the total loss of sale. Loss of sale at each station and in each hour can be defined as the difference between the total demand of bikes (number of bikes that start their trip from the station) and total supply of bikes.\n",
    "\n",
    "Total supply is comprised of number of bikes that end their trip at the station plus all the existing bikes at the station (a.k.a inventory) plus number of bikes that are added or removed from that station in that hour through some bike transfers. \n",
    "\n",
    "**Assumptions:**\n",
    "- Inventory at the beginning of first hour (in our case, hour 7) is zero.\n",
    "- At any given hour, we have access to a limited number of bikes that can be added to the stations in hope of helping reduce the imbalance without yet transferring the bikes between stations (since this analysis is during morning rush hours)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79287692",
   "metadata": {},
   "source": [
    "## Load Required Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fd07b80",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-16T00:28:00.613144Z",
     "start_time": "2022-11-16T00:28:00.583600Z"
    }
   },
   "outputs": [],
   "source": [
    "stations = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/bike_share/top_stations.csv', index_col='station')\n",
    "# run locally\n",
    "# stations = pd.read_csv('top_stations.csv', index_col='station')\n",
    "stations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d6872ea",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-16T00:28:34.586233Z",
     "start_time": "2022-11-16T00:28:34.559120Z"
    }
   },
   "outputs": [],
   "source": [
    "stations_flow = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/bike_share/stations_flow.csv')\n",
    "# run locally\n",
    "# stations_flow = pd.read_csv('stations_flow.csv')\n",
    "stations_flow['datetime'] = pd.to_datetime(stations_flow['datetime'])\n",
    "stations_flow.head(25)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a791ab30",
   "metadata": {},
   "source": [
    "The `stations_flow` data contains the prediction for the first 5 days of August 2022 and during all the hours. Our analysis is for morning rush hours, between 7 to 9 am. Also, we can run our MO model daily. For now, we only focus on the first day but at the end, we will show the full model and how it can be run daily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "190b9bce",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-16T00:30:20.677910Z",
     "start_time": "2022-11-16T00:30:20.654847Z"
    }
   },
   "outputs": [],
   "source": [
    "# Pandas will give a few SettingWithCopyWarning here when new columns are created. \n",
    "# They are false alarms. So, we suppress them.\n",
    "pd.options.mode.chained_assignment = None\n",
    "morning_flow = stations_flow[stations_flow['datetime'].dt.hour.between(7, 9)]\n",
    "morning_flow['date'] = morning_flow['datetime'].dt.date\n",
    "morning_flow['time'] = morning_flow['datetime'].dt.hour\n",
    "# For now, let's run the MO model for the first date: 08/01/2022\n",
    "flow_df = morning_flow.loc[morning_flow['date'] == datetime.date(2022, 8, 1)]\n",
    "flow_df.set_index(['station', 'time'], inplace=True)\n",
    "flow_df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75881443",
   "metadata": {},
   "source": [
    "## Problem Formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b8229be",
   "metadata": {},
   "source": [
    "Let's define our notations for the MO model. We want to run the model for every hour and every station. So, let's define the following two sets:\n",
    "\n",
    "**Sets**\n",
    "- $I\\quad$: Set of stations\n",
    "- $T\\quad$: Set of hours \n",
    "\n",
    "We also have the following information from `stations` and `flow_df` dataframes:\n",
    "\n",
    "**Parameters**\n",
    "- $e_{i,t}\\quad$: number of bikes that end their trip at station $i$ at hour $t$ (a.k.a. supply)\n",
    "- $s_{i,t}\\quad$: number of bikes that start their trip at station $i$ at hour $t$ (a.k.a. demand)\n",
    "- $c_{i}\\quad$: capacity of station $i$\n",
    "\n",
    "We know it's not easy to transfer bikes between stations during rush hours and heavy traffic. To mitigate loss due to unavailable bikes at high-demand stations, we assume there is a small reserve of bikes available at the beginning of an hour that we can allocate to the stations. We show that with $N$: \n",
    "- $N\\quad$: Number of bikes at hand that we can assign to stations at a given hour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8140631",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-15T15:51:55.077210Z",
     "start_time": "2022-11-15T15:51:55.066655Z"
    }
   },
   "outputs": [],
   "source": [
    "num_bikes = 25  # N: Number of bikes at hand that we can assign to stations at a given hour\n",
    "station_names = list(stations.index)  # set I\n",
    "time_rng = morning_flow['time'].drop_duplicates().values  # set T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68f17fc6-abf3-470d-9848-a94d98c1f7cc",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T18:23:39.820414Z",
     "start_time": "2022-11-13T18:23:39.804231Z"
    }
   },
   "outputs": [],
   "source": [
    "station_time = flow_df.index  # pair of (i,t) index\n",
    "start_forecast = flow_df.start_forecast  # s\n",
    "end_forecast = flow_df.end_forecast  # e\n",
    "capacity = stations.capacity  # c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fcf3cf8",
   "metadata": {},
   "source": [
    "## Decision Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d21f6f93",
   "metadata": {},
   "source": [
    "First of all, we like to find out how many bikes should be added to or removed from each station during each hour. We show number of bikes added with $y_{i,t}$ and those removed with $z_{i,t}$:\n",
    "- $y_{i,t}\\quad$: number of bikes to be added to station $i$ at hour $t$\n",
    "- $z_{i,t}\\quad$: number of bikes to be removed from station $i$ at hour $t$\n",
    "\n",
    "The goal of the model is to reduce the lost sale at each station per hour. This value, also depends on the value of decision variables $y_{i,t}$ and $z_{i,t}$ and as a result, is a decision variable itself. We show that with $l_{i,t}$:\n",
    "- $l_{i,t}\\quad$: lost sale at station $i$ at hour $t$\n",
    "\n",
    "**Any other variables that you can think of?**\n",
    "\n",
    "With these initial decision variables, we can start the MO model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18f9b67f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-14T18:40:14.941185Z",
     "start_time": "2022-11-14T18:40:14.311062Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "mdl = gp.Model('bike_rebalancing')\n",
    "# Variables\n",
    "\n",
    "y = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='y')\n",
    "# y = mdl.addVars(station_names, time_rng, lb=0, vtype=GRB.CONTINUOUS, name='y')  # alternatively\n",
    "z = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='z')\n",
    "l = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='l')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53a33172",
   "metadata": {},
   "source": [
    "## Constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "225d8693",
   "metadata": {},
   "source": [
    "First, we set up lowerbound and upperbound values for the number of bikes that can be added to or removed from a station in each hour.\n",
    "\n",
    "In each hour, we have $s_{i,t}$ bikes that start their trip from the station and $e_{i,t}$ bikes that end their trip at the station. If $s_{i,t} \\ge e_{i,t}$, demand is exceeding supply. In this case, we may choose to add some bikes to this station to reduce the loss (we can also add nothing). One thing to remember is that we cannot add more bikes than the station's capacity.\n",
    "\n",
    "On the other hand, if $s_{i,t} \\le e_{i,t}$, we may choose to remove some of the excess bikes from that station (we can also remove nothing). However, if $e_{i,t} \\ge s_{i,t} + c_i$, then we'll have more bikes arriving to the station than even the station capacity. In that case, we must remove some bikes to avoid overflow.\n",
    "\n",
    "The above descriptions, help us to define the bounds on $y_{i,t}$ and $z_{i,t}$.\n",
    "For adding bikes, the bounds are:\n",
    "\n",
    "\\begin{align}\n",
    "&??? &\\quad \\forall i \\in I, t \\in T\\\\\n",
    "\\end{align}\n",
    "\n",
    "For removing bikes:\n",
    "\n",
    "\\begin{align}\n",
    "&??? &\\quad \\forall i \\in I, t \\in T\\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe03af94-b72d-440c-b81d-e1de82eed7c6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-14T18:56:42.840230Z",
     "start_time": "2022-11-14T18:56:42.827230Z"
    }
   },
   "outputs": [],
   "source": [
    "# Bounds of Y\n",
    "... # ???\n",
    "\n",
    "# Bounds of z\n",
    "... # ???"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db534a31",
   "metadata": {},
   "source": [
    "Next, we set up the initial inventory (i.e. inventory at hour $t_0$) to be 0. \n",
    "\n",
    "\\begin{align}\n",
    "&??? & \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0b73d5f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T18:25:20.279863Z",
     "start_time": "2022-11-13T18:25:20.270344Z"
    }
   },
   "outputs": [],
   "source": [
    "t0 = 7\n",
    "# setting up initial inventory\n",
    "... # ???"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28688fc6",
   "metadata": {},
   "source": [
    "How do you define the inventory at a station at the beginning of an hour? Is there a limit on the amount of inventory?\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "&??? &\\\\\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5faf6bdf",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfbe6926",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1ab9a6a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "fdfeb12a",
   "metadata": {},
   "source": [
    "Next, we need to define how the loss of sale is calculated. Loss of sale is the difference between the demand and all the supply of bikes at a station.\n",
    "- Demand of bikes are all the bikes that leave a station, in any shape or form. So, what are the demands?\n",
    "- Supply of bikes are all the bikes that arrive at a station, in any shape or form. So, what are the supplies?\n",
    "\n",
    "Of course, if the supply is greater than the demand, there is no loss. So, we need to ensure that loss only considers non-negative values. This can be achieved by: \n",
    "\n",
    "\\begin{align}\n",
    "&??? &\\quad \\forall i \\in I, t \\in T \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a467dca",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:36:22.912966Z",
     "start_time": "2022-11-13T19:36:22.885648Z"
    }
   },
   "outputs": [],
   "source": [
    "# loss definition\n",
    "... # ???"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90b4f56e",
   "metadata": {},
   "source": [
    "We assumed that we have a small reserve of bikes at the beginning of each hour to allocate to stations. This limit is on the total number of bikes added to the stations. \n",
    "\n",
    "What do you think we need to do to add this constraint?\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "&??? \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61f52e65",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:36:26.781812Z",
     "start_time": "2022-11-13T19:36:26.773870Z"
    }
   },
   "outputs": [],
   "source": [
    "# limit on number of bikes added\n",
    "... # ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "904b93a6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bbd480c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a9e0d5df",
   "metadata": {},
   "source": [
    "## Objective"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cafeb700",
   "metadata": {},
   "source": [
    "The objective is to minimize total loss of sale. \n",
    "\n",
    "$$\\min \\sum_{i,t} l_{i,t}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7352d11",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:36:29.138482Z",
     "start_time": "2022-11-13T19:36:29.124154Z"
    }
   },
   "outputs": [],
   "source": [
    "mdl.setObjective(l.sum(), GRB.MINIMIZE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93782ff9",
   "metadata": {},
   "source": [
    "We can now tell Gurobi that the model is complete and it can solve the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b4669ff",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:36:30.309795Z",
     "start_time": "2022-11-13T19:36:30.214344Z"
    }
   },
   "outputs": [],
   "source": [
    "mdl.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9a1a8f4",
   "metadata": {},
   "source": [
    "## Post Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f1b8fd1",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:36:42.952139Z",
     "start_time": "2022-11-13T19:36:42.867732Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame()\n",
    "if mdl.status == GRB.Status.OPTIMAL:\n",
    "    df = flow_df.copy()\n",
    "    df = df.merge(stations[['capacity']], left_on='station', right_index=True)\n",
    "    df[['bikes_added', 'bikes_removed', 'loss_sale', 'beginning_inventory']] = 0\n",
    "    for k, v in y.items():\n",
    "        df.loc[k, 'bikes_added'] = v.x\n",
    "        df.loc[k, 'bikes_removed'] = z[k].x\n",
    "        df.loc[k, 'beginning_inventory'] = q[k].x\n",
    "        df.loc[k, 'loss_sale'] = l[k].x\n",
    "    df.reset_index(inplace=True)\n",
    "    print(f'Total Loss : {mdl.objVal}')\n",
    "    total_bikes_added = df.groupby('time')['bikes_added'].sum()\n",
    "    print(f'Total number of bikes added in each hour:\\n {total_bikes_added}')\n",
    "else:\n",
    "    print('Could not find a feasible solution!')\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94c40922",
   "metadata": {},
   "source": [
    "# Model Enhancement"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1863587a",
   "metadata": {},
   "source": [
    "Looking at the `output_df`, you may notice stations that see a transfer of bikes (addition or removal) at every hour. Even worse are those stations where some bikes are removed from them in one hour but then bikes are added to them in the next hour. We know these additions/removals consume time and money. What can we do to avoid this situation?\n",
    "\n",
    "One way to formulate this is to introduce a fixed cost for the use of the truck (or the bicycle trailer) that transfers the bikes and then add this term to the objective function. With a cost associated with the transfer, the model is incentivized to use fewer number of transfers.\n",
    "\n",
    "So first, we should calculate the number of times a transfer has occurred. Any time bikes are added to a station or removed from a station, a transfer has happened. So, we need a way to link addition of bikes and removal of bikes to a transfer. One way to achieve this is to introduce two new binary variables as follows:\n",
    "\n",
    "- $x_{i,t}\\quad$: 1 if any bike is added to station $i$ at hour $t$; 0 otherwise\n",
    "- $w_{i,t}\\quad$: 1 if any bike is removed from station $i$ at hour $t$; 0 otherwise\n",
    "\n",
    "Next, we need to establish the relationship between $y_{i,t}$ with $x_{i,t}$, and $z_{i,t}$ with $w_{i,t}$. Basically, we want to say:\n",
    "if $y_{i,t} \\ge 0$, then $x_{i,t} = 1$ and if $z_{i,t} \\ge 0$, then $w_{i,t} = 1$. \n",
    "\n",
    "We introduce the following two constraints:\n",
    "\n",
    "\\begin{align}\n",
    "&y_{i,t} \\le M x_{i,t} &\\quad \\forall i \\in I, t \\in T \\tag{11}\\\\\n",
    "&z_{i,t} \\le M w_{i,t} &\\quad \\forall i \\in I, t \\in T \\tag{12}\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $M$ is a large number.\n",
    "\n",
    "Constraint 11 ensures that if $y_{i,t}\\ge 0$, then $x_{i,t} = 1$. But by itself, this constraint cannot make $x_{i,t} =0$ if $y_{i,t} \\le 0$. The same is true with constraint 12. It ensures that $z_{i,t} \\ge 0$ make $w_{i,t} = 1$. However, it cannot force $w_{i,t} = 0$ if $z_{i,t} \\ge 0$.\n",
    "\n",
    "This can be achieved by the objective function.\n",
    "\n",
    "Total number of transfers is equal to sum of $x_{i,t}$ and $w_{i,t}$ and our goal is to minimize number of transfers. So, we add these terms to the objective function. In other words, our new objective function is:\n",
    "\n",
    "$$\\min \\sum_{i,t} (l_{i,t} + x_{i,t} + w_{i,t})$$\n",
    "\n",
    "Since minimizing total transfers is desired, the objective function tries to make both $x_{i,t}$ and $w_{i,t}$ as small as possible (or zero, in this case). Along with constraints 11 and 12, this means that for cases where $x_{i,t}$ and $w_{i,t}$ can take either 0 or 1, objective function forces them to get a value of 0. Moreover, since any extra transfer causes either $x_{i,t}$ or $w_{i,t}$ to be 1, the model is incentivized to have fewer transfers in order to minimize the objective function.\n",
    "\n",
    "Of course, you can make this model even more generic by having a coefficient for each term in the objective function (think of them as cost)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8c272b6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:37:17.713950Z",
     "start_time": "2022-11-13T19:37:17.356485Z"
    }
   },
   "outputs": [],
   "source": [
    "x = mdl.addVars(station_time, vtype=GRB.BINARY, name='x')  # 1 if y_{i,t} >= 0 \n",
    "w = mdl.addVars(station_time, vtype=GRB.BINARY, name='w')  # 1 if z_{i,t} >= 0\n",
    "\n",
    "big_m = 1000  # large number\n",
    "# relation between y and x\n",
    "mdl.addConstrs((y[i, t] <= big_m * x[i, t] for i, t in station_time), 'rel_y_x')\n",
    "# relation between z and w\n",
    "mdl.addConstrs((z[i, t] <= big_m * w[i, t] for i, t in station_time), 'rel_z_w')\n",
    "\n",
    "# new objective\n",
    "obj = l.sum() + (x.sum() + w.sum())\n",
    "mdl.setObjective(obj, GRB.MINIMIZE)\n",
    "mdl.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98a94077",
   "metadata": {},
   "source": [
    "## Post Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd6baf18",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-13T19:37:23.168043Z",
     "start_time": "2022-11-13T19:37:23.108652Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame()\n",
    "if mdl.status == GRB.Status.OPTIMAL:\n",
    "    df = flow_df.copy()\n",
    "    df = df.merge(stations[['capacity']], left_on='station', right_index=True)\n",
    "    df[['bikes_added', 'bikes_removed', 'loss_sale', 'beginning_inventory']] = 0\n",
    "    for k, v in y.items():\n",
    "        df.loc[k, 'bikes_added'] = v.x\n",
    "        df.loc[k, 'bikes_removed'] = z[k].x\n",
    "        df.loc[k, 'beginning_inventory'] = q[k].x\n",
    "        df.loc[k, 'loss_sale'] = l[k].x\n",
    "    df.reset_index(inplace=True)\n",
    "    print(f'Total Loss : {mdl.objVal}')\n",
    "    total_bikes_added = df.groupby('time')['bikes_added'].sum()\n",
    "    print(f'Total number of bikes added in each hour:\\n {total_bikes_added}')\n",
    "else:\n",
    "    print('Could not find a feasible solution!')\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cf3c78f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-10-19T18:40:17.606647Z",
     "start_time": "2022-10-19T18:40:17.596562Z"
    }
   },
   "source": [
    "# Extra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "690cd096",
   "metadata": {},
   "source": [
    "## Scenario Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50335f47",
   "metadata": {},
   "source": [
    "An important requirements in many MO problems is scenario analysis or what-if analysis. \n",
    "Generally speaking, in what-if analysis, we're interested in knowing how the solution changes under various scenarios. \n",
    "Think about our case here.\n",
    "- How does the solution change if the number of reserved bikes increase or decrease by 10%?\n",
    "- In the enhanced model, what happens if the cost of lost sale change? How about the cost of visiting a location for adding or removing the bikes?\n",
    "- What if a new station is added close to the busiest station?\n",
    "- What if we want to ensure that every station have at least 2 available bikes at the beginning of each hour?\n",
    "\n",
    "Our model here is still a simple model. But you can imagine the value that the scenario analysis can provide. It can enable you to answer many questions by creating and comparing different scenarios and evaluating their outcomes, so that you can assess their impacts on the business goals. To learn more, you can check this [multi-scenario example](https://www.gurobi.com/documentation/9.5/examples/multiscenario_py.html) from Gurobi."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57187c2f",
   "metadata": {},
   "source": [
    "## How this problem is solved in reality?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a504980e",
   "metadata": {},
   "source": [
    "After knowing how many bikes needed in each station, the bikes need to be physically moved from one station to another. \n",
    "\n",
    "During rush hours, the traffic is already heavy. So, bicycle trailers (that can usually hold 5 bicycles) are used to move bikes around.\n",
    "\n",
    "During lighter hours (mainly at night), the bikes are transferred using trucks.\n",
    "\n",
    "In either case, some bikes should be removed from stations where there are more in-flow of bikes and should be transferred to stations where there are more out-flow of bikes to balance them out. \n",
    "This problem, where trucks need to go from one station to another and either pick up bikes or deliver them, is itself another mathematical optimization problem. \n",
    "\n",
    "In this problem, we need to ensure that all the stations that have a pickup or delivery, are visited during a certain time window and the goal can be to do this with minimum number of trucks or minimum transportation cost (for example, fuel cost plus the cost for using the truck). This problem is a variation of the famous Vehicle Routing Problem (VRP) with pickup and delivery.\n",
    "To learn more, check out [this webinar](https://www.gurobi.com/resource/how-to-synchronize-complex-routing-operations-synched-vrps-with-gurobi/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e47e0051",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
